/*
   This is a very brief program in order to test the ewald code.

   In the header file 'ewald.h' there is a maximal limit for the
   reciprocal space cutoff, which (if required) has to be adjusted.
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include<iostream>

//#include "mymath.h"
//#include "ewald.h"

/*----------------------------------------------------------------------*/
double rcut  = 4.9;      /* real space cutoff. Note that this must be smaller than L/2, if minimum image convention is to be used! */
int    kmax  = 10;        /* maximal k-vector */
double Ewaldalpha=1.0;    /*Value of the Ewald Splitting parameter */


/*----------------------------------------------------------------------*/
/* number of particles: */
#define NP 8

/* sum of square of charges */
double Q2 = 8.0;

/* dielectric constant of 'surrounding' medium */
double epsilon = 1.0;   /* = vacuum */

/* Define three arrays for the particle-coordinates: */
double x[NP], y[NP], z[NP];

/* Define one array for the particle-charges: */
double Q[NP];

/* Define three arrays for the particle-forces: */
double Fx[NP], Fy[NP], Fz[NP];

/* Coulomb-energy of the box */
double E_Coulomb;

/* Length of the system */
double L = 10.0;

/* real space cutoff. Note that this must be smaller than L/2, if
   minimum image convention is to be used! */
double rcut2;   /* rcut^2 */


/* prefactor of the Coulomb potential */
double coulomb_prefactor = 1.0;

/* Berechnet das Quadrat von double x, gibt double zurueck. */
double SQR(double x) { return x*x; }

/* Wert von double x runden, gibt double zurueck. */
double dround(double x) { return floor(x+0.5); }

#define     PI     3.14159265358979323846264338328
#define   wupi     1.77245385090551602729816748334

/* constant to use when alpha is to be optimized */
#define ALPHA_OPT 0.0

 void	Ewald_init(double, int, double,
  		   double, int, double, double, double);
 double Ewald_cut(void);
 double Ewald_r_space(double *, double *, double *, double *, int,
  		      double *, double *, double *);
 double Ewald_k_space(double *, double *, double *, double *, int,
  		      double *, double *, double *);
 double Ewald_k_subset(double *, double *, double *, double *, int, int);
 double Ewald_dipol(double *, double *, double *, double *, int,
  		    double *, double *, double *);
 double Ewald_get_alpha();
 double Ewald_KPerr_k();
 double Ewald_KPerr_r();
 double Ewald_KPerr();

/*----------------------------------------------------------------------*/

int main(void)
{
  int i;

  double s1 = 0.50*L;
  double eDirecto, eReciproco;

  /* using the formula of Benson-Mackenzie on gets for the Madelung
     constant of the NaCl lattice: */
  const double madelung_NaCl = -1.747564594633182190636212035544397403481;

  /* define the particle coordinates and charges such such that this
     gives an NaCl-lattice upon periodic replication: */

  x[ 0] = 0;  y[ 0] = 0;  z[ 0] = 0;  Q[ 0] =  1.0;
  x[ 1] = s1; y[ 1] = s1; z[ 1] = 0;  Q[ 1] =  1.0;
  x[ 2] = s1; y[ 2] = 0;  z[ 2] = s1; Q[ 2] =  1.0;
  x[ 3] = 0;  y[ 3] = s1; z[ 3] = s1; Q[ 3] =  1.0;
  x[ 4] = s1; y[ 4] = 0;  z[ 4] = 0;  Q[ 4] = -1.0;
  x[ 5] = 0;  y[ 5] = s1; z[ 5] = 0;  Q[ 5] = -1.0;
  x[ 6] = 0;  y[ 6] = 0;  z[ 6] = s1; Q[ 6] = -1.0;
  x[ 7] = s1; y[ 7] = s1; z[ 7] = s1; Q[ 7] = -1.0;

  rcut2 = SQR(rcut);

  /* initialize the Ewald-algorithm */
    Ewaldalpha=ALPHA_OPT;

    Ewald_init(L,NP,Ewaldalpha,rcut,kmax,coulomb_prefactor,Q2,epsilon);

  printf("Tuned Ewald: alpha_opt=%e KPerr=%e\n", Ewald_get_alpha(), Ewald_KPerr());

  /* set forces and energy to 0 at the beginning: */
  for (i = 0; i < NP; i++) Fx[i] = Fy[i] = Fz[i] = 0.0;
  E_Coulomb = 0.0;

  /* calculate the real space contribution to the coulomb energy/forces */
  eDirecto = Ewald_r_space(x,y,z,Q,1,Fx,Fy,Fz);
  E_Coulomb += eDirecto;

  /* calculate the reciprocal space contribution to the coulomb energy/forces */
  eReciproco= Ewald_k_space(x,y,z,Q,1,Fx,Fy,Fz);
  E_Coulomb +=eReciproco;

  /*     WE DO NOT CALCULATE THE DIPOLE CONTRIBUTION !!!
	 ...but if we wished, this would be done like:
	 E_Coulomb += Ewald_dipol(x,y,z,Q,1,Fx,Fy,Fz)
  */

  /* print the results. Note that in a perfect crystal the forces on
     all particles should be identically zero! */

  for (i = 0; i < NP; i++)
    /*printf("Fx[%d]=% le\tFy[%d]=% le\tFz[%d]=% le\n",i,Fx[i],i,Fy[i],i,Fz[i])*/;
  printf("\nEnergia directa: %le\n",eDirecto);
  printf("\nEnergia reciproca: %le\n",eReciproco);
  printf("\nCoulomb-energy of the Box: %le\n",E_Coulomb);
  printf("\nCoulomb-energy of the Box: %le\n",E_Coulomb);
  printf("according to that, you predict a  Madelung constant = %14.10lf\n",L*E_Coulomb/NP);
  printf("relative error, respect real Madelung constant: %le\n",(L*E_Coulomb/NP-madelung_NaCl)/madelung_NaCl);

  /*E_Coulomb += Ewald_cut();
  printf("\n");
  printf("relative error (with correction): %le\n",
	 (L*E_Coulomb/NP-madelung_NaCl)/madelung_NaCl);

  */

 return 0;
}



//=======================================================================
//=======================================================================
//=======================================================================
//=======================================================================
//=======================================================================
//=======================================================================

/* maximal k-vector */
#define Maxkmax 200

/* maximal precision (used in cutoff) */
#define PREC 1.0e-30


/* precision of alpha */
#define ALPHA_OPT_PREC 1.0e-10


extern void   Ewald_init(double, int, double,
			 double, int, double, double, double);
extern double Ewald_cut(void);
extern double Ewald_r_space(double *, double *, double *, double *, int,
			    double *, double *, double *);
extern double Ewald_k_space(double *, double *, double *, double *, int,
			    double *, double *, double *);
extern double Ewald_k_subset(double *, double *, double *, double *, int, int);
extern double Ewald_dipol(double *, double *, double *, double *, int,
			  double *, double *, double *);
extern double Ewald_get_alpha();
extern double Ewald_KPerr_k();
extern double Ewald_KPerr_r();
extern double Ewald_KPerr();


static double  Li;       /* box length and its inverse */
static double  alpha;      /* Ewald parameter */
static double  rmax,rmax2; /* real space cutoff and its square */
static int     kmax2; /* reciprocal space cutoff and its square */
static double  prefactor;  /* prefactor of the Coulomb potential */
static double  dipfac;     /* 2*PI / ((1 + 2*epsilon) * L^3) */

/* stores the influence function: */
static double  Ghat[Maxkmax+1][Maxkmax+1][Maxkmax+1];

static double  e_re[2*Maxkmax+1][2*Maxkmax+1][2*Maxkmax+1];
static double  e_im[2*Maxkmax+1][2*Maxkmax+1][2*Maxkmax+1];

/*----------------------------------------------------------------------*/

/* internal function */

static void compute_influence_function(void);

/*----------------------------------------------------------------------*/

/* external functions */

void   Ewald_init(double length, int particlenumber, double Ewaldalpha,
		  double realspacecutoff, int reciprocalspacecutoff,
		  double coulombprefactor, double Qsquare, double epsilon);
double Ewald_cut(void);
double Ewald_r_space(double *x, double *y, double *z, double *Q, int force_flag,
		     double *Fx, double *Fy, double *Fz);
double Ewald_k_space(double *x, double *y, double *z, double *Q, int force_flag,
		     double *Fx, double *Fy, double *Fz);
double Ewald_k_subset(double *x, double *y, double *z, double *Q, int start, int end);
double Ewald_dipol(double *x, double *y, double *z, double *Q, int force_flag,
		   double *Fx, double *Fy, double *Fz);
double Ewald_KPerr_k();
double Ewald_KPerr_r();
double Ewald_KPerr();

double Ewald_get_alpha();

/*----------------------------------------------------------------------*/

static void compute_influence_function(void)
{
  /* calculates the influence function 2.0/L^2 * exp(-(PI*n/(alpha*L))^2)/n^2
     as a function of lattice vector n (NOT k=2*PI*n/L). This is stored in the
     array Ghat[Maxkmax+1][Maxkmax+1][Maxkmax+1]. For symmetrie reasons only
     one octant is actually needed. */

  int    nx,ny,nz;
  double qua,fak1,fak2;


  fak1 = 2.0/SQR(L);
  fak2 = SQR(PI/(alpha*L));

/*   fprintf(stderr,"calculating influence function for the standard Ewald method using\n" */
/* 	         "L=%lf\tkmax=%d\talpha=%lf ...",L,kmax,alpha); */

  for (nx=0; nx<=kmax; nx++)
    for (ny=0; ny<=kmax; ny++)
      for (nz=0; nz<=kmax; nz++) {
	if ((nx==0) && (ny==0) && (nz==0)) Ghat[nx][ny][nz]=0.0;
	else {
	  qua = SQR(nx) + SQR(ny) + SQR(nz);
	  Ghat[nx][ny][nz] = fak1 * exp(-fak2*qua) / qua;
	}
      }


/*   fprintf(stderr,"\n"); */
}


/*----------------------------------------------------------------------*/

double compute_KPerr_r(double alpha) {
  /* compute the real space part of the Kolafa Perram error estimate
     */
   double res;

  res = Q2*sqrt(rmax/(2.0*L*L*L)) * exp(-SQR(alpha)*rmax2) / (SQR(alpha)*rmax2);

  return res;
}

double compute_KPerr_k(double alpha) {
  /* compute the k space part of the Kolafa Perram error estimate */
  double res;

  res = Q2 * alpha
    * exp(-SQR(PI*kmax/(alpha*L)))
    / (SQR(PI) * pow(kmax, 1.5));

  return res;
}

static void compute_optimal_alpha(void) {
  /* use bisectional method to get optimal alpha value */
  double alpha_low, f_low;
  double alpha_high, f_high;
  double alpha_guess, f_guess;

  alpha_low = 0.01;
  alpha_high = 10.0;

  f_low = compute_KPerr_r(alpha_low) - compute_KPerr_k(alpha_low);
  f_high = compute_KPerr_r(alpha_high) - compute_KPerr_k(alpha_high);

  if (f_low*f_high > 0.0) {
    printf("Error: Could not init method to find optimal alpha!\n");
    exit(1);
  }

  do {
    alpha_guess = 0.5 *(alpha_low + alpha_high);
    f_guess = compute_KPerr_r(alpha_guess) - compute_KPerr_k(alpha_guess);
    if (f_low*f_guess < 0.0) alpha_high = alpha_guess;
    else alpha_low = alpha_guess;
  } while (fabs(alpha_low-alpha_high) > ALPHA_OPT_PREC);

  alpha = 0.5 *(alpha_low + alpha_high);
}


/*----------------------------------------------------------------------*/

void Ewald_init(double length, int particlenumber, double Ewaldalpha,
		double realspacecutoff, int reciprocalspacecutoff,
		double coulombprefactor, double Qsquare, double epsilon)
{

  L         = length;
  Li        = 1.0/L;
/*  NP        = particlenumber;*/
  rmax      = realspacecutoff;
  rmax2     = SQR(rmax);
  kmax      = reciprocalspacecutoff;
  kmax2     = SQR(kmax);
  prefactor = coulombprefactor;
  Q2        = Qsquare;
  dipfac    = 2.0*PI / ((1.0 + 2.0*epsilon) * L*L*L);

  /* get the alpha value */
  if (Ewaldalpha == ALPHA_OPT) compute_optimal_alpha();
  else alpha     = Ewaldalpha;

  if (kmax > Maxkmax)
    fprintf(stderr,"Subroutine 'Ewald-init' complains:\n"
	           "Desired maximal k-vector (%d) is larger\n"
	           "than what the program can cope with (%d)\n"
	           "Program terminated!\n\n",kmax,Maxkmax);

  compute_influence_function();


  /*  fprintf(stderr,"'Ewald' successfully initialized!\n"); */
}

/*----------------------------------------------------------------------*/

/* compute the cutoff error */
double Ewald_cut(void) {
	double erfc(double x);
/*   static double Cr = 0.086169460839742334501423215443; */
/*   static double Ck = 0.076533059679637957528619551795; */
  static double C = -2.837297479480619610825442578061;
  int nx, ny, nz;
  double kx, ky, kz, k2;
  double ghatsum;
  double zeta, zetacutk, zetacutr;
  double Uimgcut;
  double fak1, fak2;

  zeta = C*Li;

  // compute the sum of the influence function
  ghatsum = 0.0;
  fak1 = 2.0*PI*Li;
  fak2 = -0.25 / SQR(alpha);
  for (nx = -kmax+1; nx <= kmax; nx++) {
    kx = fak1*nx;
    for (ny = -kmax+1; ny <= kmax; ny++) {
      ky = fak1*ny;
      for (nz = -kmax+1; nz <= kmax; nz++) {
	if (nx != 0 || ny != 0 || nz != 0) {
	  if (SQR(nx) + SQR(ny) + SQR(nz) < kmax2) {
	    kz = fak1*nz;
	    k2 = SQR(kx) + SQR(ky) + SQR(kz);
	    ghatsum += 4.0*PI/k2 * exp(fak2 * k2);
	  }
	}
      }
    }
  }

  zetacutk = Li*Li*Li*ghatsum - 2.0*alpha/wupi;
  zetacutr = 2.0*PI * Li*Li*Li *
    (rmax2
     - rmax/(alpha*wupi) * exp(-SQR(alpha)*rmax2)
     - erfc(alpha*rmax) * (rmax2 + 1.0/(2.0*SQR(alpha))));

  Uimgcut = 0.5*Q2 * (zeta - zetacutk - zetacutr);

  return Uimgcut;
}

double erfc(double x)
{
	double A1 = 0.254829592;
	double A2 = -0.284496736;
	double A3 = 1.421413741;
	double A4 = -1.453152027;
	double A5 = 1.061405429;
	double P  =  0.3275911;
	double T, XSQ, TP;

	T  = 1.0/(1.0+P*x);
	XSQ = x*x;
	TP = T*(A1+T*(A2+T*(A3+T*(A4+T*A5))));
	return TP*exp(-XSQ);
}

/*----------------------------------------------------------------------*/

double Ewald_dipol(double *x, double *y, double *z, double *Q,
		   int force_flag, double *Fx, double *Fy, double *Fz)
{
  double dipx=0.0, dipy=0.0, dipz=0.0;
  double d1=prefactor*dipfac, d2=0.0;
  int    i;

  for (i=0; i<NP; i++) {
    dipx += Q[i]*x[i];
    dipy += Q[i]*y[i];
    dipz += Q[i]*z[i];
  }

  d2 += d1 * ( SQR(dipx) + SQR(dipy) + SQR(dipz) );

  if (force_flag) {
    d1   *= 2.0;
    dipx *= d1;
    dipy *= d1;
    dipz *= d1;
    for (i=0; i<NP; i++) {
      Fx[i] -= Q[i] * dipx;
      Fy[i] -= Q[i] * dipy;
      Fz[i] -= Q[i] * dipz;
    }
  }
  return d2;
}
/*----------------------------------------------------------------------*/

double Ewald_k_space(double *x, double *y, double *z, double *Q,
		     int force_flag, double *Fx, double *Fy, double *Fz)
{
  int    i;
  int    nx,ny,nz;
  double h,c,s,re,im;
  double d1,d2,d3=0.0;


  d1 = prefactor * L / (4.0 * PI);

  for (nx=-kmax; nx<=kmax; nx++)
    for (ny=-kmax; ny<=kmax; ny++)
      for (nz=-kmax; nz<=kmax; nz++)
	e_re[nx+kmax][ny+kmax][nz+kmax] = e_im[nx+kmax][ny+kmax][nz+kmax] = 0.0;

  for (nx=-kmax; nx<=kmax; nx++)
    for (ny=-kmax; ny<=kmax; ny++)
      for (nz=-kmax; nz<=kmax; nz++) {
	for (i=0; i<NP; i++)
	  if (SQR(nx) + SQR(ny) + SQR(nx) < kmax2) {
	    h = 2.0 * PI * (nx*x[i] + ny*y[i] + nz*z[i]) / L;
	    e_re[nx+kmax][ny+kmax][nz+kmax] += Q[i]*cos(h);
	    e_im[nx+kmax][ny+kmax][nz+kmax] -= Q[i]*sin(h);
	  }
	for (i=0; i<NP; i++)
	  if (SQR(nx) + SQR(ny) + SQR(nz) < kmax2) {
	    h = 2.0 * PI * (nx*x[i] + ny*y[i] + nz*z[i]) / L;

	    c = cos(h);
	    s = sin(h);

	    d2 = Ghat[abs(nx)][abs(ny)][abs(nz)];
	    re = (e_re[nx+kmax][ny+kmax][nz+kmax]*c
		- e_im[nx+kmax][ny+kmax][nz+kmax]*s) * d2;

	    if (force_flag) {
	      im = (e_re[nx+kmax][ny+kmax][nz+kmax]*s
                  + e_im[nx+kmax][ny+kmax][nz+kmax]*c) * d2;

	      d2 = im * Q[i];

	      Fx[i] += nx * d2;
	      Fy[i] += ny * d2;
	      Fz[i] += nz * d2;
	    }
	    d3 += Q[i] * re;
	  }
      }
  d3 *= d1;

  /* self energy correction: */
  d3 -= prefactor * Q2 * alpha / wupi;

  return d3;
}

/*----------------------------------------------------------------------*/

double Ewald_k_subset(double *x, double *y, double *z, double *Q, int start, int end)
{
  /* calculates the contribution to the Ewald k-space energy resulting from
     interactions of charges with numbers in {start,start+1,...,end}.
     Note that 'start' and 'end' are INCLUDED in the list! */
  int    i;
  int    nx,ny,nz;
  int    absnx, absny, absnz;
  int    nxpkmax,nypkmax,nzpkmax;
  double k,d1,d2,d3=0.0;
  double Q2_subset=0.0;


  d1 = prefactor * L / (4.0 * PI);
  d2 = 2.0 * PI / L;

  for (i = start; i <= end; Q2_subset += SQR(Q[i++]));

  for (nx = -kmax; nx <= kmax; nx++)
    for (ny = -kmax; ny <= kmax; ny++)
      for (nz = -kmax; nz <= kmax; nz++)
	e_re[nx+kmax][ny+kmax][nz+kmax] = e_im[nx+kmax][ny+kmax][nz+kmax] = 0.0;

  for (nx = -kmax; nx <= kmax; nx++) {
    absnx = abs(nx);
    nxpkmax = nx + kmax;
    for (ny = -kmax; ny <= kmax; ny++) {
      absny = abs(ny);
      nypkmax = ny + kmax;
      for (nz = -kmax; nz <= kmax; nz++) {
	absnz = abs(nz);
	nzpkmax = nz + kmax;
	for (i = start; i <= end; i++)
	  if (SQR(nx) + SQR(ny) + SQR(nx) < kmax2) {
	    k = d2 * (nx*x[i] + ny*y[i] + nz*z[i]);
	    e_re[nxpkmax][nypkmax][nzpkmax] += Q[i] * cos(k);
	    e_im[nxpkmax][nypkmax][nzpkmax] -= Q[i] * sin(k);
	  }
	for (i = start; i <= end; i++)
	  if (SQR(nx) + SQR(ny) + SQR(nx) < kmax2) {
	    k = d2 * (nx*x[i] + ny*y[i] + nz*z[i]);
	    d3 += Q[i] * Ghat[absnx][absny][absnz] *
	          ( e_re[nxpkmax][nypkmax][nzpkmax] * cos(k)
		  - e_im[nxpkmax][nypkmax][nzpkmax] * sin(k) );
	  }
      }
    }
  }
  d3 *= d1;

  /* self energy correction: */
  d3 -= prefactor * Q2_subset * alpha / wupi;

  return d3;
}

/*----------------------------------------------------------------------*/

double Ewald_r_space(double *x, double *y, double *z, double *Q,
		     int force_flag, double *Fx, double *Fy, double *Fz)
{
  int    t1,t2;
  double dx,dy,dz,r,r2,ar,ar2;
  double d1,d2,d3=0.0,d4;
  double erfc(double x);

  for (t1 = 0; t1 < NP-1; t1++)   /* Quick and Dirty N^2 loop */
    for (t2 = t1 + 1; t2<NP; t2++) {
      dx = x[t1] - x[t2]; dx -= dround(dx*Li)*L;
      dy = y[t1] - y[t2]; dy -= dround(dy*Li)*L;
      dz = z[t1] - z[t2]; dz -= dround(dz*Li)*L;
      r2 = SQR(dx) + SQR(dy) + SQR(dz);
	if (r2 <= rmax2) {
	  ar2 = SQR(ar = alpha * (r = sqrt(r2)));
	  d3 += ( d2 = (d1 = Q[t1] * Q[t2] * prefactor) * erfc(ar) / r );

	  if (force_flag) {
	    d4  = ( d2 + d1*2.0*alpha*exp(-ar2)/wupi) / r2;

	    Fx[t1] += d4*dx;
	    Fy[t1] += d4*dy;
	    Fz[t1] += d4*dz;
	    Fx[t2] -= d4*dx;
	    Fy[t2] -= d4*dy;
	    Fz[t2] -= d4*dz;
	  }
	}
    }
  return d3;
}

/*----------------------------------------------------------------------*/

double Ewald_KPerr_r() {
  return compute_KPerr_r(alpha);
}

double Ewald_KPerr_k() {
  return compute_KPerr_k(alpha);
}

double Ewald_KPerr() {
  double KPerr_r, KPerr_k;

  /* compute the Kolafa-Perram error */
  KPerr_r = compute_KPerr_r(alpha);
  KPerr_k = compute_KPerr_k(alpha);
  return KPerr_r + KPerr_k;
}

double Ewald_get_alpha() { return alpha; }



















